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^ ■ ABSTRACT 

X5 : . 

1^ ^ High-Q optical resonances in photonic microcavities are investigated numerically using a time-harmonic finitc- 
pLn ' element method. 
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I r 1. INTRODUCTION 

C/2 ■ 

O Optical microcavities allow to confine light to small volumes. High resonance Q-factors can be attained using the 
high reflectivity of multi-layer Fabry-Perot resonators, of total internal reflection and/or of photonic bandgap 
materials. '^'^l Jn this contribution we revisit a microcavity design proposed originally by Notomi, Kuramochi and 
Taniyama.^ Notomi et al have shown that very high Q-factors can be reached with a size-modulated ID stack 
^ design. In this case confinement to the cavity is obtained by total internal reflection in two dimensions, and by 
^ ' a photonic bandgap in the third dimension. Advantages of this design are its compactness and simplicity. 

For an efficient design of microcavities and other integrated photonic components 3D simulations of Maxwell's 
' equations are needed. We have developed finite-element method (FEM) based solvers for the Maxwell eigenvalue 
and for the Maxwell scattering problems. The method is based on higher order vectorial elements, adaptive 
unstructured grids, and on a rigorous treatment of transparent boundaries. We perform a numerical analysis of 
the microcavity setup. Results on resonance wavelengths and Q-factors are generated using a resonance solver. 
The obtained values are confirmed by simulations of transmission spectra of light incident to the microcavities. 
Dependence of Q-factors and resonance wavelength on structural parameters are investigated. A focus of this 
contribution is on the numerical convergence of the obtained results. 

^ : 

(Sj ; 2. INVESTIGATED SETUP 

■ The investigated setup consists of a finite array of dielectric blocks, as depicted in Figure [TJ The first and the 
I . last block of the array are merged to waveguides with same square shape as the blocks, but infinite in the third 

• • ' dimension (z-direction). The centers of the blocks are arranged periodically with a period of a = 430 nm. Width 

■ s-nd height of the waveguide and blocks (in x- and y-direction) are Wx = 3a and Wy = 0.5a. The length of 
^ ■ 2Ni + 1 inner blocks is modulated symmetrically around a central block, Wz^n = Wzfl + {\n\/Niy{Wz,Ni — W^2,o), 

5^ . with Wzfi = 0.45a, Wz,n, = 0.5345a. The length of 2No outer blocks is constant Wz,n = Wz,Ni (n > Ni), 
where the distance from the center block is denoted with n (in number of blocks) . The setup of Notomi et al 
corresponds to Ni = 13. The blocks consist of dielectric with refractive index n = 3.46, the surrounding is filled 
with vacuum (n = 1.0). Conceptually similar devices have also been investigated experimentally.!^ 
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3. NUMERICAL RESULTS 



3.1 Band diagram of the ID periodic array 

We start the numerical investigation by computing the band diagram of the ID periodic array of dielectric blocks 
in vacuum. The unit cell of for the computation is periodic in z-direction with a period a, transparent boundary 
conditions (perfectly matched layers, PML) apply in x- and y-directions, the width of the blocks in z-direction is 
Wz = 0.5345 a. Computation of modes with even symmetry is selected by reducing the symmetric computational 
domain to the half space a: > and applying corresponding boundary conditions at the symmetry plane. Figure[2] 
(left) shows a part of the band structure of the ID photonic crystal.'^ A band gap appears for dimensionless 
frequencies between Q « 0.246 and Cj « 0.329. Dimensionless frequency Cj is defined as a) = (a;a)/(27rco), with 
time-harmonic frequency uj and vacuum speed of light cq. Dimensionless magnitude of Bloch vector k is defined 
as k — fca/(27r), with Bloch vector k. For a ~ 430 nm the bandgap corresponds to a wavelength range between 
approximately 1307nm and 1748 nm. In Figure [5] (left) only modes below the light cone are displayed. Modes 
above the light cone have complex eigenfrequencies, which corresponds to leakage to the transversal directions.'^ 

Figured] (right) shows numerical convergence of the simulated lowest Bloch-mode at fc = 0.5: Finite-element 
simulations have been performed for the same physical parameters and for different numerical parameters with 
increasing numerical accuracy. In this case the polynomial degree of the finite-element ansatzfunctions,® p, has 
been varied between p — 1 and p = 7 for a fixed spatial mesh discretizing the geometry of the unit cell. The 
computed value of Hj for p = 1 , ujp=7 ~ 0.24644, is taken as quasi-exact value, and the relative error of the 
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Figure 1. Top: Geometry of a size-modulated ID stack microcavity with a central block of width Wzfi, 2 x A^i = 2 x 13 
size-modulated blocks of widths Wz,n and 2 x No = 2x5 blocks of constant width Wz,Ni- Blocks and incoming/outgoing 
waveguides consist of dielectric material, the structure is surrounded by air, waveguides extend to infinity. The structure 
is mirror-symmetric at the planes x = 0, j/ = 0, therefore the computational domain can be reduced to one quarter of the 
whole structure. Bottom: Field distribution on resonance (pseudo-color representations), from top to bottom: intensity 
in a a;-z-cross-section, intensity in a j/-z-cross-section, '^(Ex) in a x-z-cross-section, K(ii'i^) in a y-2-cross-section, log(_E^) 
in a a;-2-cross-section, log(_E^) in a y-z-cross-section. 
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Figure 2. Left: Part of the band structure of the ID photonic crystal. The hght cone is indicated by a grey hne. Right: 
Convergence of the relative error of the computed eigenfrequency of the lowest Bloch mode at fc = 0.5. 



computed values of w for p = 1 ... 6, Acjp — \ojp — u)p^7\/ujp=r is plotted as function of the number of unknowns 
of the discrete problem. Computation times on a standard workstation range between below a second and few 
minutes for the plotted values. High accuracy with relative errors better than 0.1% is obtained for p > 3. The 
solvers used for the band diagram simulations and for the further simulations throughout this contribution are 
part of the FEM program package JCMsuit^ which is developed by Zuse Institute Berlin and JCMwave. 

3.2 Transmission of waveguide modes through a finite ID periodic array 

For simulating transmission through a finite array of blocks / a finite ID photonic crystal at specific wavelength 
Aq, we first compute the fundamental propagation mode of the waveguide at Aq using a FEM propagation mode- 
solver. The obtained mode field is applied as input data to one of the boundaries of the 3D computational domain 
(left boundary in Fig. [T]), such that the mode propagates in direction of the center of the waveguide. We then 
compute the scattered light field in the setup corresponding to this excitation using higher-order finite-elements. 
In post-processes we extract energy fluxes through interfaces and field distributions in several cross-sections from 
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Figure 3. Transmission spectrum of incident waveguide modes through a finite ID periodic array of 15 blocks of length 
Wz = 0.5345 a. Transmission is suppressed in the frequency range between u ~ 0.246 and to ~ 0.329. 



the 3D light field distribution. Transmission T is defined as ratio of the energy fiux of the outgoing light field 
at the back side of the array to the energy flux of the incoming waveguide mode through the waveguide cross 
section. 

Figure [3] shows a transmission spectrum for 800 incident waveguide modes with vacuum wavelengths between 
Aq = 1075 nm (w = 0.4) and Aq = 2150 nm (w = 0.2). In perfect agreement with the band structure simulations 
transmission is greatly suppressed in the frequency range between Co w 0.246 and Cj w 0.329 which corresponds 
to a bandgap of the ID periodic structure ( c/ Figure [5]) . 

3.3 Direct simulation of high-Q cavity resonances 

We use an eigenmode solver for directly simulating the resonance properties (resonance wavelength and Q-factor) 
of size-modulated ID stack cavities: Given the geometrical setup as described in Section [2l one computes an 
electric held distribution E and a complex eigenfrequency w which satisfy Maxwell's time-harmonic wave equation 
without sources, 

V X ^"^V xE = uj^sE , 

electric permittivity and magnetic permeability are denoted by e and fi, respectively. Transparent boundary 
conditions (PML) take into account the specific geometry of the problem where waveguides are modelled to 
extend to infinity in the exterior domain. When the eigenmode (E, uj) is computed, the respective Q-factor is 
deduced from the real and imaginary parts of the complex eigenfrequency, Q — 3?(cj)/(— 23(cj)), the resonance 
wavelength Aq is given by Ao = 27rco/K(aj). Figure [T] shows visualizations of a typical field distribtion obtained 
with the resonance solver. 

Direct simulation of a resonance requires a single computation only. In contrast, deducing the resonance 
properties from a transmission or reflection spectrum generated with a time-harmonic solver requires several 
computations at various wavelengths. Especially for high-Q resonances, where the choice of wavelengths for 
a transmission scan is not clear a-priori, direct computation of resonances simplifies the simulation task and 
reduces computational effort. Time-domain solvers for simulating 3D high-Q resonances typically suffer from 
very long computation times and slow convergence rates. 

We have performed simulations for different cavity setups. Figure H] shows the dependence of Q-factor and 
resonance wavelength on the number of modulated blocks, Ni, of the ID stack cavity setup (see Section [2|). 
For these simulations, the number of outer blocks is fixed to No = 4. With increasing number of modulated 
blocks, the resonance wavelength is pulled from the band-edge (at a wavelength of 1748 nm) more and more 
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Figure 4. Dependence of resonance wavelength Ao (left) and Q-factor (right) on number of modulated blocks, Ni. For all 
simulations. No — 4. Results are displayed for finite-element degree p = 3 and p — 4. 



0.02 
0.018 
0.016 
0.014 

'Fo.oia 

o 0.01 
<l 0.008 
0.006 
0.004 
0.002 


Figure 5. 
No. 



10 



4 6 

No 



' O " o o o 



10 



10 : 



u 

I 

a 



10 : 



10 



No 



10 



Dependence of resonance wavelength Ao (left) and Q-factor (right) on number of outer, unmodulated blocks, 



inside the band-gap. At the same time, the Q-factor increases nearly exponentially for the investigated range 
of A^i = 1 . . . 15. Figure 0] shows results for different numerical resolutions, i.e., for finite-element degree p = 3 
and p = 4. As can be seen from these results, for higher numbers of modulated blocks, i.e., for higher Q-factors, 
influence of numerical resolution on the results gets more significant, up to a level of relative deviations of the 
order of ten percent. 

Figure [S] shows the dependence of Q-factor and resonance wavelength on the number of outer, unmodulated 
blocks, No, of the ID stack cavity setup (see Section [5]). For these simulations, the number of inner blocks is 
fixed to Ni = 13 which corresponds to the setting of Notomi et afP For these simulations, fourth-order finite- 
elements {p = 4) have been used. With increasing number of unmodulated blocks, the resonance wavelength is 
changed only on a sub-nanometer scale. The Q-factor increases exponentially for A"o < 6 and reaches a plateau 
of Q w 3 X 10^ for No > 7. We expect that radiation losses to the sourrounding air limit high-Q performance in 
this case.l^ 

For evaluating numerical errors of high-Q resonance simulations using the resonance mode solver we again 
show a numerical convergence study: The the geometry setting is fixed to to Ni = 13, No — 5. Finite-element 
polynomial degree p is varied between p — 1 and p = 5, and Q-factor, resonance wavelength and numerical effort 
are recorded. Figure [5] (left) shows how the relative error of Q-factor and resonance wavelength converge with 
number of unknowns N, where the solution at highest resolution (p — 5) has been taken as reference result 
(Age = 1694.12 nm, Qqe — 1.04 x 10^, Nqe = 978,570). Computation times for the displayed results are 1 sec, 
12 sec, 2min, 9min, respectively (on a workstation with 8 CPU cores). This demonstrates that resonances with 
Q-factors of about 1 million can be computed at accuracies of about 0.1 % (resonance wavelength) and 10% (Q) 
within few minutes. 

In principle, a wrong numerical realization of transparent boundary conditions can introduce errors which 
will not show up in a convergence study as shown in Figure [6] (left). We realize transparent boundary conditions 
with the PML method, i.e., by a coordinate transform to complex space, by FEM discretization with higher-order 
elements (same order p as in the interior domain), and by truncation at a variable distance from the boundary 
to the interior computational domain. For a fixed setting of p = 4 we have further investigated dependence of 
the results on varied number of PML segments. Figure [5] (right) shows that convergence of the results with PML 
discretization parameter is very fast. This indicates that the main contribution to the numerical errors as shown 
in Figure[6] (left) are caused by discretization of the electromagnetic field in the interior domain. Typical settings 
for the results displayed in Figure H] (left) and Figure SI [5] are about seven PML segments. The setting of further 
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Figure 6. Convergence of numerical results for high-Q (10^) resonance simulation using a resonance solver: Dependence 
of relative errors of Q-factor and resonance wavelength, AQ — \Qn — Qqe\/Qqe and A A = |Ajv — Xqe\/Xqe, on number of 
unknowns of the FEM problem, N. Geometry setting: A*'; = 13, A'^o — 5. Left: Finite-element degree p varied between 
p = 1 and p = 4. Right: Number of segments discretizing the solution in the PML region varied between 1 and 7. 



PML parameters is done adaptivelyP More detailed studies should also consider influence of these parameters 
on numerical accuracy. 

We note that an evaluation of the significant differences between the numerical results of Notomi et aW and 
our numerical results requires further investigations. Notomi et al have reported results from a finite-difference 
time-domain method yielding significantly higher Q-factors (by about two orders of magnitude) for the same 
physical setting. 

3.4 Light scattering response of a high-Q cavity 

An alternative possibility of computing the resonance properties of a cavity is to investigate its scattering response 
to incident light fields. Using the same approach as in Section we perform FEM light scattering simulations 
where the incident light fields are waveguide modes at fixed frequencies/ wavelengths. From the near field solutions 
we deduce the energy stored in the central part of the cavity, and the energy flux transmitted through the device. 

For these simulations we again investigate the modulated stack cavity with Ni = 13, No — 5. As numerical 
discretization parameter we choose a polynomial degree of the finite-element ansatzfunctions of p = 4. Figure [7] 
(left) shows the (normalized) energy E stored in the cavity as a function of the wavelength of the incident 
waveguide mode. For automatic determination of the wavelengths of the incoming waveguide modes, Ai„ we use 
a self-adaptive bisection algorithm. In this example we used a search range of AA ~ 3 nm and 14 bisections, 
resulting in a final wavelength-resolution of 5\ ~ 10""* nm. With a typical simulation time of about lOmin for a 
single- wavelength simulation, this results in a total computation time for the frequency scan of about 8 hrs. 

We fit a Lorentzian distribution to the data points and deduce a resonance wavelength of Aq ~ 1694.7345 nm 
and a Q-factor of Q w 1.14 x 10^ [Xi^/FWHM). Both are in very good agreement with the results from the 
resonance mode solver within the ranges of numerical error (see Section l3.3p . as expected from the convergence 
study. 

We have also recorded the total transmission through the device by integrating the energy fluxes over the 
incoming and outgoing facets of the device. Figure [7] (right) shows the transmission peak corresponding to the 
resonance. Here the total transmission reaches a maximum of T w 0.31. 
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Figure 7. Left: Energy stored in the cavity, E, as function of the wavelength of the incident waveguide mode, Af„. Circles: 
simulation results at distinct frequencies, line: Lorentz-fit. Right: Transmission through the cavity, T, as function of Aj„. 
Circles: simulation results at distinct frequencies, line: spline-fit. 



4. CONCLUSION 

Size- modulated stack mieroeavities have been investigated numerically using time-harmonic finite-element solvers. 
The results have been validated by convergence analysis. Well converged results for cavities with Q-factors larger 
than 10^ have been obtained. Dependencies of resonance wavelength and quality factor on geometrical parameters 
have been investigated. 
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